% This code simulates the dynamics of the Ramsey solution in a representative-agent economy.
% We consider a shock on G


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%        clearing data         %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
warning('off','all');
isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0;
isOctave && warning('off', 'Octave:array-as-logical');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%   Loading steady-state data    %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load ../toDynare % load the file produced by the Julia code


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%    Defining model variables    %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Persistence of the shock
rho_u = 0.1;

% Normalization of the standard deviation for keeping unchanged the NPV of G 
u0 = .01*(1-rho_u/(solution.R))/solution.G; 

% Importing some variables from toDynare.mat and simplifying notation
tauk   = eco.tauk;
alpha  = eco.alpha(1);
beta   = eco.beta;
delta  = eco.delta(1);
kappa  = 0;
chi    = eco.chi;
phi    = eco.phi;
G      = solution.G;
coefB  = 0;
GsY    = solution.G/solution.Y;
chi    = 1;
BsY    = 0;
dif    = 1;
r      = 1/beta - 1;
FK     = r;
KsL    = ((FK+delta)/alpha)^(1/(alpha - 1));
FL     = (1-alpha)*KsL^alpha;
w      = FL;
Ltot   =  (chi*w)^phi;
K      = KsL*Ltot;
Y      = K^alpha*Ltot^(1-alpha);
I      = delta*K;
Ctot   = Y*(1-GsY) -I; 
G      = GsY*Y;
A      = (Ctot +G - w*Ltot)/r; % check A=K
x      = Ctot - (1/chi)*Ltot^(1+1/phi)/(1+1/phi);

% utility function
if eco.sigma==1
    util = @(c) log(c);
else
    util = @(c) (c^(1-eco.sigma)-1)/(1-eco.sigma);
end

% Welfare  
W = util(x);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%   Opening the Dynare file    %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

file = 'code_dynare_RA_rules.mod';
fid  = fopen(file, 'w');
formatSpec = '%25.20f';


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%     Variables and params     %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['var\n'] ; fprintf(fid, str);

str = ['u r w FK FL K Ltot Ctot Y I It\n'] ; fprintf(fid, str);
str = ['ut Kt Ct wt rt Yt Lt G x W;\n'] ; fprintf(fid, str);

str = ['varexo eps; \n\n'] ; fprintf(fid, str);
str = ['parameters\n'] ; fprintf(fid, str);
str = ['beta alpha phi delta gamma chi rho_u Gss;\n\n'] ; fprintf(fid, str);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%        Initialization        %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['alpha','   = ',num2str(alpha,formatSpec),';\n']; fprintf(fid, str);
str = ['beta','    = ',num2str(beta,formatSpec),';\n']; fprintf(fid, str);
str = ['phi','     = ',num2str(eco.phi,formatSpec),';\n']; fprintf(fid, str);
str = ['delta','   = ',num2str(eco.delta,formatSpec),';\n']; fprintf(fid, str);
str = ['gamma','   = ',num2str(eco.sigma,formatSpec),';\n']; fprintf(fid, str);
str = ['rho_u','   = ',num2str(rho_u,formatSpec),';\n']; fprintf(fid, str); %HERE
str = ['Gss','     = ',num2str(G,formatSpec),';\n']; fprintf(fid, str);
str = ['chi','     = ',num2str(chi,formatSpec),';\n']; fprintf(fid, str);

equa = 1;  % Numbering equations
tol  = 0;  % threshold for transition, to avoid to consider very small transitions
str  = ['\n\n'] ; fprintf(fid, str);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%       Model definition       %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['model;\n\n']; fprintf(fid, str);

% Defining the variable x
str = ['x = Ctot - (1/chi)*Ltot^(1+1/phi)/(1+1/phi); %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

% Euler equation
str = ['(x^-gamma) =  beta*(1+r(+1))*(x(+1)^-gamma); %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

% Public spending shock
str = ['G = Gss*(1+u); %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

% Capital dynamics
str = ['I = K - (1-delta)*K(-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

% Resource constraint
str = ['G  + Ctot + I = Y;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

% Production side
str = ['FL = (1 - alpha)*(K(-1)/Ltot)^alpha; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['FK = alpha*(K(-1)/Ltot)^(alpha-1)-delta; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['u = rho_u*u(-1) + eps;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Y = K(-1)^alpha*Ltot^(1-alpha);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

% Wage
str = ['Ltot =  (chi*w)^phi; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

% Taxes
str = ['r = FK;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['w = FL;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

% Welfare
if eco.sigma==1
    str = ['W = log(x);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
else
    str = ['W =  (x',num2str(h),'^(1-gamma)-1)/(1-gamma) ;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
end;


%%%%%%%% Log-deviation variables
% There are denoted with 't': for variable x, xt = (x -x_ss)/x_ss, where x_ss is the steady state value of x.

str = ['ut = 100*u;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Kt = 100*(K/',num2str(K,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Ct = 100*(Ctot/',num2str(Ctot,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['wt = 100*(w/',num2str(w,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['rt = 100*(r - (1/beta - 1));  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Lt = 100*(Ltot/',num2str(Ltot,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Yt = 100*(Y/',num2str(Y,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['It = 100*(I/',num2str(delta*K,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['end;\n\n'] ; fprintf(fid, str);

str = ['%% end of the model part\n\n'] ; fprintf(fid, str);

%%%%%%%% end of the model part


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%     Steady-state model       %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['initval;\n\n'] ; fprintf(fid, str);
    
str = ['Ltot  = ', num2str(Ltot,formatSpec),';\n'] ; fprintf(fid, str);
str = ['r     = 1/beta - 1;\n'] ; fprintf(fid, str);
str = ['w     = ', num2str(w,formatSpec),';\n'] ; fprintf(fid, str);
str = ['FK    = ', num2str(FK,formatSpec),';\n'] ; fprintf(fid, str);
str = ['FL    = ', num2str(FL,formatSpec),';\n'] ; fprintf(fid, str);
str = ['K     = ', num2str(K,formatSpec),';\n'] ; fprintf(fid, str);
str = ['x     = ', num2str(x,formatSpec),';\n'] ; fprintf(fid, str);
str = ['u     = ', num2str(0),';\n'] ; fprintf(fid, str);
str = ['Ctot  = ', num2str(Ctot,formatSpec),';\n'] ; fprintf(fid, str);
str = ['Y     = ', num2str(Y,formatSpec),';\n'] ; fprintf(fid, str);
str = ['I     = delta*K;\n'] ; fprintf(fid, str);
str = ['ut    = 0;\n'] ; fprintf(fid, str);
str = ['Kt    = 0;\n'] ; fprintf(fid, str);
str = ['Ct    = 0;\n'] ; fprintf(fid, str);
str = ['wt    = 0;\n'] ; fprintf(fid, str);
str = ['rt    = 0;\n'] ; fprintf(fid, str);
str = ['Lt    = 0;\n'] ; fprintf(fid, str);
str = ['Yt    = 0;\n'] ; fprintf(fid, str);
str = ['It    = 0;\n'] ; fprintf(fid, str);

str = ['G     = ', num2str(G,formatSpec),';\n'] ; fprintf(fid, str);
str = ['W     = ', num2str(W,formatSpec),';\n'] ; fprintf(fid, str);

str = ['end;\n\n'] ; fprintf(fid, str);

%%%%%%%% end of the steady state part


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%   Steady-state computation    %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['resid;\n\n'] ; fprintf(fid, str);
str = ['steady(nocheck);\n\n'] ; fprintf(fid, str);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%   Defining aggregate shock   %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['shocks;\n var eps; stderr ',num2str(u0,formatSpec),';\n end;\n'] ; fprintf(fid, str);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%     Running  stoch_simul     %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['stoch_simul (order=1,irf=500) ut Yt Ct It Lt Kt W;'] ; fprintf(fid, str);
isOctave && fflush(fid);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%       Launching Dynare       %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

command = ['dynare ',file, ' noclearall'];
eval(command);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%         Saving data          %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

name = ['fig',num2str(rho_u),'FB.mat'];
save (name,'ut_eps','Yt_eps','It_eps','Lt_eps','Ct_eps','Kt_eps','W','W_eps','Ctot','Ltot')
